I began creating a figure by taking the fGSEA results and coloring them by pathway types, e.g. glucose vs lipid metabolism, etc. in Illustrator. Preliminary results look really good, so I want to generate these graphs through R here.
I’ll run DESeq2 → fGSEA on all timepoints together, then every individual timepoints, then take Hallmark and KEGG data for WT SPF vs GF data to generate figures. If this works, I can do the same thing for other comparisons as well, e.g. SPF WT vs KO.
For this run, I grouped pathways into categories, some broad and others specific:
Since I’ve run this DESeq2 → fGSEA pipeline many times before in the same manner on this dataset, I will skip many of the visualizations/graphs for quality control I’ve generated in the past for this.
Libraries imported.
knitr::opts_chunk$set(echo = TRUE)
klippy::klippy(position = c("top", "right"))
reqpkg <- c("feather", "DESeq2", "dplyr", "magrittr", "genefilter", "AnnotationHub", "org.Mm.eg.db", "GO.db", "vsn", "pheatmap", "biomaRt", "curl", "RColorBrewer", "viridis", "fgsea", "tidyverse", "DT", "ggplot2", "ggpubr")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "feather"
## [1] '0.3.5'
## [1] "DESeq2"
## [1] '1.24.0'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "magrittr"
## [1] '1.5'
## [1] "genefilter"
## [1] '1.66.0'
## [1] "AnnotationHub"
## [1] '2.16.1'
## [1] "org.Mm.eg.db"
## [1] '3.8.2'
## [1] "GO.db"
## [1] '3.8.2'
## [1] "vsn"
## [1] '3.52.0'
## [1] "pheatmap"
## [1] '1.0.12'
## [1] "biomaRt"
## [1] '2.40.5'
## [1] "curl"
## [1] '4.3'
## [1] "RColorBrewer"
## [1] '1.1.2'
## [1] "viridis"
## [1] '0.5.1'
## [1] "fgsea"
## [1] '1.10.1'
## [1] "tidyverse"
## [1] '1.3.0'
## [1] "DT"
## [1] '0.13'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
theme_set(theme_pubr())
select <- dplyr::select
mart <- useDataset("mmusculus_gene_ensembl", mart=useMart("ensembl"))
mouseProteinCodingGenes <- getBM(attributes=c("ensembl_gene_id","external_gene_name","description"), filters='biotype', values=c('protein_coding'), mart=mart)
bm <- getBM(attributes=c("ensembl_gene_id", "hsapiens_homolog_associated_gene_name"), mart=mart) %>%
distinct() %>%
as_tibble() %>%
na_if("") %>%
na.omit()
pathway_folder <- '../pathways/'
f_res <- list()
anno <- list()
plot_res <- list()
dict <- list()
dict$immune <- c(
"IMMUNE|INFLAMMATORY|TNFA|ALLOGRAFT|JAK_STAT|TGF_BETA|IL2_STAT5|INTERFERON_GAMMA|REACTIVE_OXYGEN_SPECIES|CYTOKINE|CHEMOKINE_SIGNALING|FC_EPSILON_RI|HEMATOPOIETIC_CELL|IMMUNODEFICIENCY|GRAFT_VERSUS_HOST|NOD_LIKE_RECEPTOR|TOLL_LIKE_RECEPTOR|FC_GAMMA_R|LEUKOCYTE|NATURAL_KILLER_CELL|T_CELL|COAGULATION|JNK_CASCADE|INTERLEUKIN|B_CELL|CELLULAR_DEFENSE_RESPONSE|LYMPHOCYTE|JUN_KINASE|CHEMOTAXIS|MACROPHAGE|MAST_CELL|MEGAKARYOCYTE|LYMPHOID_PROGENITOR_CELL|EXTRAVASATION|EOSINOPHIL|DENDRITIC_CELL|NEUTROPHIL|TUMOR_NECROSIS_FACTOR|RESPONSE_TO_BACTERIUM|DEFENSE_RESPONSE|COMPLEMENT_ACTIVATION|ANTIGEN|INTERFERON_ALPHA|T_HELPER|PLASMINOGEN"
)
dict$cancer <- c(
"CANCER|CARCINOMA|LEUKEMIA|REGULATION_OF_CELL_MORPHOGENESIS|VASCULAR_ENDOTHELIAL_GROWTH_FACTOR|CELL_GROWTH|MORPHOGENESIS|DIFFERENTIATION|PROLIFERATION|HOMEOSTASIS_OF_NUMBER_OF_CELLS|CELL_FATE_COMMITMENT|TISSUE_HOMEOSTASIS|RESPONSE_TO_GROWTH_FACTOR|UV_RESPONSE_DN|MELANOMA|DNA_REPAIR|DNA_DAMAGE_RESPONSE|DNA_DOUBLE_STRAND_BREAK|REGULATION_OF_CELL_DEATH"
)
dict$ox_phos <- c(
"OXIDATIVE_PHOSPHORYLATION|MITOCHONDRIAL_PROTON_TRANSPORTING_ATP_SYNTHASE_COMPLEX_ASSEMBLY|ELECTRON_TRANSPORT|RESPIRATORY_CHAIN|ENERGY_COUPLED_PROTON_TRANSPORT|ATP_BIOSYNTHETIC|NUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|NUCLEOSIDE_BISPHOSPHATE_METABOLIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|CELLULAR_RESPIRATION|MITOCHONDRIAL_MEMBRANE_POTENTIAL|CITRATE_CYCLE|TCA_CYCLE|TRICARBOXYLIC_ACID|POSITIVE_REGULATION_OF_ATPASE"
)
dict$glucose <- c(
"GLUCOSE|GLYCOLYSIS|GLUCONEOGENESIS|GLYCOGENESIS|GLYCOGENOLYSIS"
)
dict$carb <- c(
"GLUCOSE|GLYCOLYSIS|GLUCONEOGENESIS|GLYOXYLATE_AND_DICARBOXYLATE_METABOLISM|PYRUVATE_METABOLISM|CARBOHYDRATE_DERIVATIVE_BIOSYNTHETIC_PROCESS|MONOSACCHARIDE_BIOSYNTHETIC_PROCESS|CARBOHYDRATE_DERIVATIVE_METABOLIC_PROCESS|GLYCOSYL_COMPOUND_BIOSYNTHETIC|GLYCOSIDE_METABOLIC|GLYCOSYL_COMPOUND_CATABOLIC|RIBOSE_PHOSPHATE_BIOSYNTHETIC|GLYCOGENESIS|GLYCOGENOLYSIS|PENTOSE_PHOSPHATE_PATHWAY|CARBOHYDRATE_METABOLIC"
)
dict$lipid <- c(
"FATTY_ACID|CHOLESTEROL_HOMEOSTASIS|BILE_ACID|BIOSYNTHESIS_OF_UNSATURATED_FATTY_ACIDS|TERPENOID_BACKBONE_BIOSYNTHESIS|STEROID_BIOSYNTHESIS|ADIPOGENESIS|PHOSPHOLIPID_METABOLIC_PROCESS|LIPID_METABOLIC|LIPID_BIOSYNTHETIC|PHOSPHATIDYLCHOLINE_BIOSYNTHETIC|LIPID_HOMEOSTASIS|STEROL_METABOLIC|ACYLGLYCEROL_METABOLIC|STEROL_HOMEOSTASIS|CHOLINE_CATABOLIC|ISOPRENOID_METABOLIC|STEROID_METABOLIC|LIPOPROTEIN_METABOLIC|STEROID_BIOSYNTHETIC|LIPID_CATABOLIC|TRIGLYCERIDE_METABOLIC|LONG_CHAIN_FATTY_ACYL_COA_BIOSYNTHETIC|LIPID_OXIDATION|FATTY_ACYL_COA_BIOSYNTHETIC|CHOLESTEROL_BIOSYNTHETIC|STEROL_BIOSYNTHETIC|STEROID_HORMONE_BIOSYNTHESIS|DIACYLGLYCEROL_METABOLIC|GLYCOSPHINGOLIPID_BIOSYNTHETIC|TERPENOID_METABOLIC|LIPID_STORAGE"
)
dict$bile <- c(
"BILE_ACID"
)
dict$peroxisome <- c(
"PEROXISOME|PPAR_SIGNALING_PATHWAY"
)
dict$protein <- c(
"LYSINE_DEGRADATION|SELENOAMINO_ACID_METABOLISM|GLUTATHIONE_METABOLISM|TYROSINE_METABOLISM|ARGININE_AND_PROLINE_METABOLISM|CYSTEINE_AND_METHIONINE_METABOLISM|TRYPTOPHAN_METABOLISM|GLYCINE_SERINE_AND_THREONINE_METABOLISM|ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM|BETA_ALANINE_METABOLISM|PROTEASOME|VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION|REGULATION_OF_PROTEIN_CATABOLIC_PROCESS|PROTEOLYSIS|PROTEIN_CATABOLIC|PEPTIDE_METABOLIC|GLUTAMINE_FAMILY_AMINO_ACID_METABOLIC_PROCESS|TYROSINE_CATABOLIC|AMINO_ACID_CATABOLIC|AMINE_CATABOLIC|AMINO_ACID_FAMILY_METABOLIC|AMINO_ACID_METABOLIC|LYSINE_CATABOLIC|LIPOPROTEIN_METABOLIC|AMINO_ACID_BIOSYNTHETIC|GLUTAMATE_CATABOLIC|HOMOCYSTEINE_METABOLIC|L_PHENYLALANINE|AMINE_METABOLIC|PEPTIDE_BIOSYNTHETIC|PHENYLALANINE_METABOLISM|CARNITINE_BIOSYNTHETIC|ARGININE_BIOSYNTHETIC|GLUTATHIONE_METABOLIC|GLUTATHIONE_DERIVATIVE_BIOSYNTHETIC|TRYPTOPHAN_METABOLIC|PROTEIN_OXIDATION"
)
dict$nucleic <- c(
"PYRIMIDINE_METABOLISM|RNA_METABOLIC|RNA_CATABOLIC|NUCLEOSIDE_CATABOLIC|PYRIMIDINE_CONTAINING_COMPOUND_METABOLIC|PYRIDINE_CONTAINING_COMPOUND_METABOLIC|TRNA_THREONYLCARBAMOYLADENOSINE_METABOLIC|NUCLEOBASE_BIOSYNTHETIC|PYRIDINE_CONTAINING_COMPOUND_BIOSYNTHETIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_BIOSYNTHETIC|MRNA_CATABOLIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_METABOLIC|NUCLEOSIDE_PHOSPHATE_METABOLIC|PURINE_CONTAINING_COMPOUND_METABOLIC|NUCLEOSIDE_PHOSPHATE_BIOSYNTHETIC|PURINE_CONTAINING_COMPOUND_BIOSYNTHETIC|RIBONUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|CGMP_BIOSYNTHETIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_CATABOLIC"
)
dict$scfa <- c(
"BUTANOATE_METABOLISM|PROPANOATE_METABOLISM|BUTYRATE|ACETATE|GLUTAMATE|PROPRIONATE|ACETOACETIC_ACID"
)
dict$metabolism <- c(
"OXIDATIVE_PHOSPHORYLATION|GLUCOSE|GLUCOSE|FATTY_ACID|BILE_ACID|PEROXISOME|LYSINE_DEGRADATION|PYRIMIDINE_METABOLISM|BUTANOATE_METABOLISM|XENOBIOTIC_METABOLISM|MITOCHONDRIAL_PROTON_TRANSPORTING_ATP_SYNTHASE_COMPLEX_ASSEMBLY|GLYCOLYSIS|GLYCOLYSIS|CHOLESTEROL_HOMEOSTASIS|PPAR_SIGNALING_PATHWAY|SELENOAMINO_ACID_METABOLISM|RNA_METABOLIC|PROPANOATE_METABOLISM|LIMONENE_AND_PINENE_DEGRADATION|ELECTRON_TRANSPORT|GLUCONEOGENESIS|GLUCONEOGENESIS|BILE_ACID|GLUTATHIONE_METABOLISM|RNA_CATABOLIC|BUTYRATE|XENOBIOTIC_METABOLIC|RESPIRATORY_CHAIN|GLYCOGENESIS|GLYOXYLATE_AND_DICARBOXYLATE_METABOLISM|BIOSYNTHESIS_OF_UNSATURATED_FATTY_ACIDS|TYROSINE_METABOLISM|NUCLEOSIDE_CATABOLIC|ACETATE|METABOLIC|ENERGY_COUPLED_PROTON_TRANSPORT|GLYCOGENOLYSIS|PYRUVATE_METABOLISM|TERPENOID_BACKBONE_BIOSYNTHESIS|ARGININE_AND_PROLINE_METABOLISM|PYRIMIDINE_CONTAINING_COMPOUND_METABOLIC|GLUTAMATE|METABOLISM|ATP_BIOSYNTHETIC|CARBOHYDRATE_DERIVATIVE_BIOSYNTHETIC_PROCESS|STEROID_BIOSYNTHESIS|CYSTEINE_AND_METHIONINE_METABOLISM|PYRIDINE_CONTAINING_COMPOUND_METABOLIC|PROPRIONATE|CATABOLIC|NUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|MONOSACCHARIDE_BIOSYNTHETIC_PROCESS|ADIPOGENESIS|TRYPTOPHAN_METABOLISM|TRNA_THREONYLCARBAMOYLADENOSINE_METABOLIC|ACETOACETIC_ACID|CATABOLISM|NUCLEOSIDE_BISPHOSPHATE_METABOLIC|CARBOHYDRATE_DERIVATIVE_METABOLIC_PROCESS|PHOSPHOLIPID_METABOLIC_PROCESS|GLYCINE_SERINE_AND_THREONINE_METABOLISM|NUCLEOBASE_BIOSYNTHETIC|ANABOLIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|GLYCOSYL_COMPOUND_BIOSYNTHETIC|LIPID_METABOLIC|ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM|PYRIDINE_CONTAINING_COMPOUND_BIOSYNTHETIC|ANABOLISM|CELLULAR_RESPIRATION|GLYCOSIDE_METABOLIC|LIPID_BIOSYNTHETIC|BETA_ALANINE_METABOLISM|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_BIOSYNTHETIC|BIOSYNTHETIC|MITOCHONDRIAL_MEMBRANE_POTENTIAL|GLYCOSYL_COMPOUND_CATABOLIC|PHOSPHATIDYLCHOLINE_BIOSYNTHETIC|PROTEASOME|MRNA_CATABOLIC|CELLULAR_RESPIRATION|CITRATE_CYCLE|RIBOSE_PHOSPHATE_BIOSYNTHETIC|LIPID_HOMEOSTASIS|VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION|AEROBIC_RESPIRATION|TCA_CYCLE|GLYCOGENESIS|STEROL_METABOLIC|REGULATION_OF_PROTEIN_CATABOLIC_PROCESS|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_METABOLIC|RESPIRATORY_GASEOUS_EXCHANGE|TRICARBOXYLIC_ACID|GLYCOGENOLYSIS|ACYLGLYCEROL_METABOLIC|PROTEOLYSIS|NUCLEOSIDE_PHOSPHATE_METABOLIC|RESPONSE_TO_GLUCAGON|POSITIVE_REGULATION_OF_ATPASE|PENTOSE_PHOSPHATE_PATHWAY|STEROL_HOMEOSTASIS|PROTEIN_CATABOLIC|CARBOHYDRATE_METABOLIC|CHOLINE_CATABOLIC|PEPTIDE_METABOLIC|PURINE_CONTAINING_COMPOUND_METABOLIC|ISOPRENOID_METABOLIC|GLUTAMINE_FAMILY_AMINO_ACID_METABOLIC_PROCESS|NUCLEOSIDE_PHOSPHATE_BIOSYNTHETIC|STEROID_METABOLIC|TYROSINE_CATABOLIC|PURINE_CONTAINING_COMPOUND_BIOSYNTHETIC|LIPOPROTEIN_METABOLIC|AMINO_ACID_CATABOLIC|RIBONUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|STEROID_BIOSYNTHETIC|AMINE_CATABOLIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|LIPID_CATABOLIC|AMINO_ACID_FAMILY_METABOLIC|CGMP_BIOSYNTHETIC|TRIGLYCERIDE_METABOLIC|AMINO_ACID_METABOLIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_CATABOLIC|LONG_CHAIN_FATTY_ACYL_COA_BIOSYNTHETIC|LYSINE_CATABOLIC|LIPID_OXIDATION|LIPOPROTEIN_METABOLIC|FATTY_ACYL_COA_BIOSYNTHETIC|AMINO_ACID_BIOSYNTHETIC|CHOLESTEROL_BIOSYNTHETIC|GLUTAMATE_CATABOLIC|STEROL_BIOSYNTHETIC|HOMOCYSTEINE_METABOLIC|STEROID_HORMONE_BIOSYNTHESIS|L_PHENYLALANINE|DIACYLGLYCEROL_METABOLIC|AMINE_METABOLIC|GLYCOSPHINGOLIPID_BIOSYNTHETIC|PEPTIDE_BIOSYNTHETIC|TERPENOID_METABOLIC|PHENYLALANINE_METABOLISM|LIPID_STORAGE|CARNITINE_BIOSYNTHETIC|ARGININE_BIOSYNTHETIC|GLUTATHIONE_METABOLIC|GLUTATHIONE_DERIVATIVE_BIOSYNTHETIC|TRYPTOPHAN_METABOLIC|PROTEIN_OXIDATION"
)
filter_fgsea <- function(selection, data, col1="#800000", col2="#D6D6CE") {
on <- data %>% select(pathway) %>% filter(grepl(selection, pathway)) %>% pull()
off <- data$pathway %>% setdiff(on)
on <- data.frame(p=on, col=rep(col1, length(on)), in_category=rep(TRUE, length(on)))
off <- data.frame(p=off, col=rep(col2, length(off)), in_category=rep(FALSE, length(off)))
return(rbind(on, off))
}
create_palette <- function(x) setNames(as.character(x$p), x$col)
create_plot <- function(data, pal, title) {
p <- ggplot(data, aes(x = reorder(pathway, NES), y = NES, fill = pathway)) +
geom_col() +
coord_flip() +
labs(title = title) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5),
axis.title = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()) +
scale_fill_manual(values = names(pal), limits = pal)
return(p)
}
DESeq_run <- function(data) {
condition_list <- data.frame(row.names=colnames(df),
Genotype=rep(c(rep("WT",3), rep("KO",3)),2),
Condition=c(rep("SPF",6),rep("GF",6)))
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=data,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
dds <- estimateSizeFactors(dds)
dds$Group <- factor(paste0(dds$Genotype, dds$Condition))
dds$Group <- relevel(dds$Group, "WTSPF")
design(dds) <- ~ Group
dds <- DESeq(dds)
res <- results(dds,contrast=c("Group", "WTSPF", "WTGF"), tidy = TRUE)
return(res)
}
This is analysis using raw featurecount data.
filePath = '../../../'
df <- read.csv(paste0(filePath, 'KF_RNASeq_counts_filtered.txt')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
set_rownames(.$Geneid) %>% select(-c(Geneid,X))
condition_list <- data.frame(row.names=colnames(df),
Genotype=rep(c(rep("WT",18), rep("KO",18)),2),
Condition=c(rep("SPF",36),rep("GF",36)))
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
dds <- estimateSizeFactors(dds)
dds$Group <- factor(paste0(dds$Genotype, dds$Condition))
dds$Group <- relevel(dds$Group, "WTSPF")
design(dds) <- ~ Group
dds <- DESeq(dds)
WT_res <- results(dds,contrast=c("Group", "WTSPF", "WTGF"), tidy = TRUE)
res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
select(hsapiens_homolog_associated_gene_name, stat) %>%
na.omit() %>%
distinct() %>%
group_by(hsapiens_homolog_associated_gene_name) %>%
summarize(stat=mean(stat))
ranks <- deframe(res)
f_res$time_all$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05) %>%
arrange(desc(NES))
f_res$time_all$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
f_res$time_all$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
title <- ""
anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_all$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_all$hm <- lapply(pal, function(x) create_plot(f_res$time_all$hm, x, title))
anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_all$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_all$KEGG <- lapply(pal, function(x) create_plot(f_res$time_all$KEGG, x, title))
anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_all$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_all$GO <- lapply(pal, function(x) create_plot(f_res$time_all$GO, x, title))
filePath = '../../../data/R/by_time/'
df <- read_feather(paste0(filePath, 'time_1_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
condition_list <- data.frame(row.names=colnames(df),
Genotype=rep(c(rep("WT",3), rep("KO",3)),2),
Condition=c(rep("SPF",6),rep("GF",6)))
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
dds <- estimateSizeFactors(dds)
dds$Group <- factor(paste0(dds$Genotype, dds$Condition))
dds$Group <- relevel(dds$Group, "WTSPF")
design(dds) <- ~ Group
dds <- DESeq(dds)
WT_res <- results(dds,contrast=c("Group", "WTSPF", "WTGF"), tidy = TRUE)
res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
select(hsapiens_homolog_associated_gene_name, stat) %>%
na.omit() %>%
distinct() %>%
group_by(hsapiens_homolog_associated_gene_name) %>%
summarize(stat=mean(stat))
ranks <- deframe(res)
f_res$time_1$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05) %>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.84% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_1$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.84% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_1$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.84% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 1"
anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_1$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_1$hm <- lapply(pal, function(x) create_plot(f_res$time_1$hm, x, title))
anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_1$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_1$KEGG <- lapply(pal, function(x) create_plot(f_res$time_1$KEGG, x, title))
anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_1$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_1$GO <- lapply(pal, function(x) create_plot(f_res$time_1$GO, x, title))
df <- read_feather(paste0(filePath, 'time_2_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)
res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
select(hsapiens_homolog_associated_gene_name, stat) %>%
na.omit() %>%
distinct() %>%
group_by(hsapiens_homolog_associated_gene_name) %>%
summarize(stat=mean(stat))
ranks <- deframe(res)
f_res$time_2$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05) %>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.09% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_2$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.09% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_2$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.09% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 2"
anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_2$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_2$hm <- lapply(pal, function(x) create_plot(f_res$time_2$hm, x, title))
anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_2$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_2$KEGG <- lapply(pal, function(x) create_plot(f_res$time_2$KEGG, x, title))
anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_2$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_2$GO <- lapply(pal, function(x) create_plot(f_res$time_2$GO, x, title))
df <- read_feather(paste0(filePath, 'time_3_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)
res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
select(hsapiens_homolog_associated_gene_name, stat) %>%
na.omit() %>%
distinct() %>%
group_by(hsapiens_homolog_associated_gene_name) %>%
summarize(stat=mean(stat))
ranks <- deframe(res)
f_res$time_3$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05) %>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_3$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_3$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 3"
anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_3$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_3$hm <- lapply(pal, function(x) create_plot(f_res$time_3$hm, x, title))
anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_3$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_3$KEGG <- lapply(pal, function(x) create_plot(f_res$time_3$KEGG, x, title))
anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_3$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_3$GO <- lapply(pal, function(x) create_plot(f_res$time_3$GO, x, title))
df <- read_feather(paste0(filePath, 'time_4_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)
res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
select(hsapiens_homolog_associated_gene_name, stat) %>%
na.omit() %>%
distinct() %>%
group_by(hsapiens_homolog_associated_gene_name) %>%
summarize(stat=mean(stat))
ranks <- deframe(res)
f_res$time_4$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05) %>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.79% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_4$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.79% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_4$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.79% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 4"
anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_4$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_4$hm <- lapply(pal, function(x) create_plot(f_res$time_4$hm, x, title))
anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_4$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_4$KEGG <- lapply(pal, function(x) create_plot(f_res$time_4$KEGG, x, title))
anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_4$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_4$GO <- lapply(pal, function(x) create_plot(f_res$time_4$GO, x, title))
df <- read_feather(paste0(filePath, 'time_5_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)
res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
select(hsapiens_homolog_associated_gene_name, stat) %>%
na.omit() %>%
distinct() %>%
group_by(hsapiens_homolog_associated_gene_name) %>%
summarize(stat=mean(stat))
ranks <- deframe(res)
f_res$time_5$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05) %>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.81% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_5$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.81% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_5$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.81% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 5"
anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_5$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_5$hm <- lapply(pal, function(x) create_plot(f_res$time_5$hm, x, title))
anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_5$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_5$KEGG <- lapply(pal, function(x) create_plot(f_res$time_5$KEGG, x, title))
anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_5$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_5$GO <- lapply(pal, function(x) create_plot(f_res$time_5$GO, x, title))
df <- read_feather(paste0(filePath, 'time_6_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)
res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
select(hsapiens_homolog_associated_gene_name, stat) %>%
na.omit() %>%
distinct() %>%
group_by(hsapiens_homolog_associated_gene_name) %>%
summarize(stat=mean(stat))
ranks <- deframe(res)
f_res$time_6$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05) %>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_6$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_6$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
as_tibble() %>%
filter(padj < 0.05)%>%
arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 6"
anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_6$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_6$hm <- lapply(pal, function(x) create_plot(f_res$time_6$hm, x, title))
anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_6$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_6$KEGG <- lapply(pal, function(x) create_plot(f_res$time_6$KEGG, x, title))
anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_6$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_6$GO <- lapply(pal, function(x) create_plot(f_res$time_6$GO, x, title))
f_res$time_all$hm %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_all$KEGG %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_all$GO %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_1$hm %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_1$KEGG %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_1$GO %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_2$hm %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_2$KEGG %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_2$GO %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_3$hm %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_3$KEGG %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_3$GO %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_4$hm %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_4$KEGG %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_4$GO %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_5$hm %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_5$KEGG %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_5$GO %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_6$hm %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_6$KEGG %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
f_res$time_6$GO %>%
select(-leadingEdge, -ES, -nMoreExtreme) %>%
DT::datatable(filter = "bottom")
ggarrange(plot_res$time_all$hm$immune, plot_res$time_all$hm$cancer, plot_res$time_all$hm$metabolism, plot_res$time_all$hm$ox_phos, plot_res$time_all$hm$glucose, plot_res$time_all$hm$carb, plot_res$time_all$hm$lipid, plot_res$time_all$hm$bile, plot_res$time_all$hm$peroxisome, plot_res$time_all$hm$protein, plot_res$time_all$hm$nucleic, plot_res$time_all$hm$scfa, ncol = 6, nrow = 2,
labels = c("immune", "cancer", "all metabolism", "oxidative phosphorylation", "glucose", "carb metabolism", "lipid metabolism", "bile", "peroxisome", "protein metabolism", "nucleic acid metabolism", "SCFA metabolism"),
font.label = list(size=10, face="plain"), hjust = 0.5, label.x = 0.5) %>%
annotate_figure(top = text_grob("Hallmark - all timepoints", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_all$KEGG$immune, plot_res$time_all$KEGG$cancer, plot_res$time_all$KEGG$metabolism, plot_res$time_all$KEGG$ox_phos, plot_res$time_all$KEGG$glucose, plot_res$time_all$KEGG$carb, plot_res$time_all$KEGG$lipid, plot_res$time_all$KEGG$bile, plot_res$time_all$KEGG$peroxisome, plot_res$time_all$KEGG$protein, plot_res$time_all$KEGG$nucleic, plot_res$time_all$KEGG$scfa, ncol = 6, nrow = 2,
labels = c("immune", "cancer", "all metabolism", "oxidative phosphorylation", "glucose", "carb metabolism", "lipid metabolism", "bile", "peroxisome", "protein metabolism", "nucleic acid metabolism", "SCFA metabolism"),
font.label = list(size=10, face="plain"), hjust = 0.5, label.x = 0.5) %>%
annotate_figure(top = text_grob("KEGG - all timepoints", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_all$GO$immune, plot_res$time_all$GO$cancer, plot_res$time_all$GO$metabolism, plot_res$time_all$GO$ox_phos, plot_res$time_all$GO$glucose, plot_res$time_all$GO$carb, plot_res$time_all$GO$lipid, plot_res$time_all$GO$bile, plot_res$time_all$GO$peroxisome, plot_res$time_all$GO$protein, plot_res$time_all$GO$nucleic, plot_res$time_all$GO$scfa, ncol = 6, nrow = 2,
labels = c("immune", "cancer", "all metabolism", "oxidative phosphorylation", "glucose", "carb metabolism", "lipid metabolism", "bile", "peroxisome", "protein metabolism", "nucleic acid metabolism", "SCFA metabolism"),
font.label = list(size=10, face="plain"), hjust = 0.5, label.x = 0.5) %>%
annotate_figure(top = text_grob("GO BP - all timepoints", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$immune, plot_res$time_2$hm$immune, plot_res$time_3$hm$immune, plot_res$time_4$hm$immune, plot_res$time_5$hm$immune, plot_res$time_6$hm$immune) %>%
annotate_figure(top = text_grob("Hallmark – immune", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$immune, plot_res$time_2$KEGG$immune, plot_res$time_3$KEGG$immune, plot_res$time_4$KEGG$immune, plot_res$time_5$KEGG$immune, plot_res$time_6$KEGG$immune) %>%
annotate_figure(top = text_grob("KEGG – immune", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$immune, plot_res$time_2$GO$immune, plot_res$time_3$GO$immune, plot_res$time_4$GO$immune, plot_res$time_5$GO$immune, plot_res$time_6$GO$immune) %>%
annotate_figure(top = text_grob("GO BP – immune", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$cancer, plot_res$time_2$hm$cancer, plot_res$time_3$hm$cancer, plot_res$time_4$hm$cancer, plot_res$time_5$hm$cancer, plot_res$time_6$hm$cancer) %>%
annotate_figure(top = text_grob("Hallmark – cancer", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$cancer, plot_res$time_2$KEGG$cancer, plot_res$time_3$KEGG$cancer, plot_res$time_4$KEGG$cancer, plot_res$time_5$KEGG$cancer, plot_res$time_6$KEGG$cancer) %>%
annotate_figure(top = text_grob("KEGG – cancer", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$cancer, plot_res$time_2$GO$cancer, plot_res$time_3$GO$cancer, plot_res$time_4$GO$cancer, plot_res$time_5$GO$cancer, plot_res$time_6$GO$cancer) %>%
annotate_figure(top = text_grob("GO BP – cancer", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$metabolism, plot_res$time_2$hm$metabolism, plot_res$time_3$hm$metabolism, plot_res$time_4$hm$metabolism, plot_res$time_5$hm$metabolism, plot_res$time_6$hm$metabolism) %>%
annotate_figure(top = text_grob("Hallmark – all metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$metabolism, plot_res$time_2$KEGG$metabolism, plot_res$time_3$KEGG$metabolism, plot_res$time_4$KEGG$metabolism, plot_res$time_5$KEGG$metabolism, plot_res$time_6$KEGG$metabolism) %>%
annotate_figure(top = text_grob("KEGG – all metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$metabolism, plot_res$time_2$GO$metabolism, plot_res$time_3$GO$metabolism, plot_res$time_4$GO$metabolism, plot_res$time_5$GO$metabolism, plot_res$time_6$GO$metabolism) %>%
annotate_figure(top = text_grob("GO BP – all metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$ox_phos, plot_res$time_2$hm$ox_phos, plot_res$time_3$hm$ox_phos, plot_res$time_4$hm$ox_phos, plot_res$time_5$hm$ox_phos, plot_res$time_6$hm$ox_phos) %>%
annotate_figure(top = text_grob("Hallmark – oxidative phosphorylation", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$ox_phos, plot_res$time_2$KEGG$ox_phos, plot_res$time_3$KEGG$ox_phos, plot_res$time_4$KEGG$ox_phos, plot_res$time_5$KEGG$ox_phos, plot_res$time_6$KEGG$ox_phos) %>%
annotate_figure(top = text_grob("KEGG – oxidative phosphorylation", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$ox_phos, plot_res$time_2$GO$ox_phos, plot_res$time_3$GO$ox_phos, plot_res$time_4$GO$ox_phos, plot_res$time_5$GO$ox_phos, plot_res$time_6$GO$ox_phos) %>%
annotate_figure(top = text_grob("GO BP – oxidative phosphorylation", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$glucose, plot_res$time_2$hm$glucose, plot_res$time_3$hm$glucose, plot_res$time_4$hm$glucose, plot_res$time_5$hm$glucose, plot_res$time_6$hm$glucose) %>%
annotate_figure(top = text_grob("Hallmark – glucose metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$glucose, plot_res$time_2$KEGG$glucose, plot_res$time_3$KEGG$glucose, plot_res$time_4$KEGG$glucose, plot_res$time_5$KEGG$glucose, plot_res$time_6$KEGG$glucose) %>%
annotate_figure(top = text_grob("KEGG – glucose metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$glucose, plot_res$time_2$GO$glucose, plot_res$time_3$GO$glucose, plot_res$time_4$GO$glucose, plot_res$time_5$GO$glucose, plot_res$time_6$GO$glucose) %>%
annotate_figure(top = text_grob("GO BP – glucose metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$carb, plot_res$time_2$hm$carb, plot_res$time_3$hm$carb, plot_res$time_4$hm$carb, plot_res$time_5$hm$carb, plot_res$time_6$hm$carb) %>%
annotate_figure(top = text_grob("Hallmark – carbohydrate metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$carb, plot_res$time_2$KEGG$carb, plot_res$time_3$KEGG$carb, plot_res$time_4$KEGG$carb, plot_res$time_5$KEGG$carb, plot_res$time_6$KEGG$carb) %>%
annotate_figure(top = text_grob("KEGG – carbohydrate metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$carb, plot_res$time_2$GO$carb, plot_res$time_3$GO$carb, plot_res$time_4$GO$carb, plot_res$time_5$GO$carb, plot_res$time_6$GO$carb) %>%
annotate_figure(top = text_grob("GO BP – carbohydrate metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$lipid, plot_res$time_2$hm$lipid, plot_res$time_3$hm$lipid, plot_res$time_4$hm$lipid, plot_res$time_5$hm$lipid, plot_res$time_6$hm$lipid) %>%
annotate_figure(top = text_grob("Hallmark – lipid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$lipid, plot_res$time_2$KEGG$lipid, plot_res$time_3$KEGG$lipid, plot_res$time_4$KEGG$lipid, plot_res$time_5$KEGG$lipid, plot_res$time_6$KEGG$lipid) %>%
annotate_figure(top = text_grob("KEGG – lipid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$lipid, plot_res$time_2$GO$lipid, plot_res$time_3$GO$lipid, plot_res$time_4$GO$lipid, plot_res$time_5$GO$lipid, plot_res$time_6$GO$lipid) %>%
annotate_figure(top = text_grob("GO BP – lipid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$bile, plot_res$time_2$hm$bile, plot_res$time_3$hm$bile, plot_res$time_4$hm$bile, plot_res$time_5$hm$bile, plot_res$time_6$hm$bile) %>%
annotate_figure(top = text_grob("Hallmark – bile acid", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$bile, plot_res$time_2$KEGG$bile, plot_res$time_3$KEGG$bile, plot_res$time_4$KEGG$bile, plot_res$time_5$KEGG$bile, plot_res$time_6$KEGG$bile) %>%
annotate_figure(top = text_grob("KEGG – bile acid", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$bile, plot_res$time_2$GO$bile, plot_res$time_3$GO$bile, plot_res$time_4$GO$bile, plot_res$time_5$GO$bile, plot_res$time_6$GO$bile) %>%
annotate_figure(top = text_grob("GO BP – bile acid", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$peroxisome, plot_res$time_2$hm$peroxisome, plot_res$time_3$hm$peroxisome, plot_res$time_4$hm$peroxisome, plot_res$time_5$hm$peroxisome, plot_res$time_6$hm$peroxisome) %>%
annotate_figure(top = text_grob("Hallmark – peroxisome", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$peroxisome, plot_res$time_2$KEGG$peroxisome, plot_res$time_3$KEGG$peroxisome, plot_res$time_4$KEGG$peroxisome, plot_res$time_5$KEGG$peroxisome, plot_res$time_6$KEGG$peroxisome) %>%
annotate_figure(top = text_grob("KEGG – peroxisome", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$peroxisome, plot_res$time_2$GO$peroxisome, plot_res$time_3$GO$peroxisome, plot_res$time_4$GO$peroxisome, plot_res$time_5$GO$peroxisome, plot_res$time_6$GO$peroxisome) %>%
annotate_figure(top = text_grob("GO BP – peroxisome", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$protein, plot_res$time_2$hm$protein, plot_res$time_3$hm$protein, plot_res$time_4$hm$protein, plot_res$time_5$hm$protein, plot_res$time_6$hm$protein) %>%
annotate_figure(top = text_grob("Hallmark – protein/amino acid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$protein, plot_res$time_2$KEGG$protein, plot_res$time_3$KEGG$protein, plot_res$time_4$KEGG$protein, plot_res$time_5$KEGG$protein, plot_res$time_6$KEGG$protein) %>%
annotate_figure(top = text_grob("KEGG – protein/amino acid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$protein, plot_res$time_2$GO$protein, plot_res$time_3$GO$protein, plot_res$time_4$GO$protein, plot_res$time_5$GO$protein, plot_res$time_6$GO$protein) %>%
annotate_figure(top = text_grob("GO BP – protein/amino acid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$nucleic, plot_res$time_2$hm$nucleic, plot_res$time_3$hm$nucleic, plot_res$time_4$hm$nucleic, plot_res$time_5$hm$nucleic, plot_res$time_6$hm$nucleic) %>%
annotate_figure(top = text_grob("Hallmark – nucleic acid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$nucleic, plot_res$time_2$KEGG$nucleic, plot_res$time_3$KEGG$nucleic, plot_res$time_4$KEGG$nucleic, plot_res$time_5$KEGG$nucleic, plot_res$time_6$KEGG$nucleic) %>%
annotate_figure(top = text_grob("KEGG – nucleic acid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$nucleic, plot_res$time_2$GO$nucleic, plot_res$time_3$GO$nucleic, plot_res$time_4$GO$nucleic, plot_res$time_5$GO$nucleic, plot_res$time_6$GO$nucleic) %>%
annotate_figure(top = text_grob("GO BP – nucleic acid metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$hm$scfa, plot_res$time_2$hm$scfa, plot_res$time_3$hm$scfa, plot_res$time_4$hm$scfa, plot_res$time_5$hm$scfa, plot_res$time_6$hm$scfa) %>%
annotate_figure(top = text_grob("Hallmark – SCFA metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$KEGG$scfa, plot_res$time_2$KEGG$scfa, plot_res$time_3$KEGG$scfa, plot_res$time_4$KEGG$scfa, plot_res$time_5$KEGG$scfa, plot_res$time_6$KEGG$scfa) %>%
annotate_figure(top = text_grob("KEGG – SCFA metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
ggarrange(plot_res$time_1$GO$scfa, plot_res$time_2$GO$scfa, plot_res$time_3$GO$scfa, plot_res$time_4$GO$scfa, plot_res$time_5$GO$scfa, plot_res$time_6$GO$scfa) %>%
annotate_figure(top = text_grob("GO BP – SCFA metabolism", face = "bold", size = 14),
left = text_grob("pathways", rot = 90),
bottom = text_grob("Normalized Enrichment Score"))
I generated the items in the categories by taking the lists of results and filtering key words and phrases specific to the category that can be used to pull all results containing said phrase. For example, for “glucose” one phrase was glycolysis.
To save time, after every analysis I checked only the unique results to sort into categories.
There is greater immune-related function in SPF mice which can be attributed to the essential role of microbes in immune development. Surprisingly, there is also increased cancer-related pathways in SPF mice as well, and GF mice have enrichment in tumor-preventing pathways, which suggests a role between microbes and oncogene expression. This may be related to the inflammation model of cancer, where any type of chronic immune response drives oncogene function, or this may highlight the relationship and balance between healthy immune functioning and oncogenesis.
Furthermore, there is broadly increased metabolic processes in GF mice accross all types of macromolecules. A minor note on the results, while it may look like protein metabolism is the most enriched for metabolism in GF mice, this may be more a reflection of the many amino acid related pathways that might all actually be centrally connected. Nonetheless, the extreme increase in host metabolism can reflect several possibilities: the host is compensating for reduced absorptive and metabolic capabilities in the gut due to the lack of microbes, or without microbes key regulatory components for metabolic homeostasis is dysregulated, and the mice are forced to rely on systemic metabolic increase to fulfill basic energy requirements, or core physiologic defects due to impaired developement in the GF condition, such as reduced mitochondrial function or brown fat is leading to upregulation of metabolic processes. With regard to the latter, in light of the decreased host immunity and cellular growth processes, an alternative hypothesis could be that the mice simply have more energy to spend in the absence of immune or heavy mitotic demands.
Digging deeper into the increased metabolic functioning of GF mice, it seems they rely more on lipid and protein sources for metabolism rather than glucose, or generally carbohydrates. Perhaps this reflects decreased carbohydrate digestion and absorption without bacteria, as other groups have shown bacteria are important for metabolite production from many starches and complex sugars.
In any case, oxidative phosphorylation specifically is the most enriched pathway in GF mice in nearly every case, thus, cellular respiration is, without a doubt, higher in GF mice. I believe this might be a culmination effect of increased metabolism accross the board — i.e. metabolism of all macromolecules are feeding into the electron tranport chain in the end to generate useable energy. It is important to note that these samples are liver samples, so I can’t comment on oxidative phosphorylation versus uncoupled electron transport energy generation in thermogenesis, but it is clear that mitochondria are playing a key role in regulating energy production of GF mice.
In the future, I would want to include another category for mitochondria, since many mitochondrial pathways showed up and are interesting, such as apoptotic or mitochondrial RNA expression pathways. Also, I want to take a look at the leading edge genes for important and recurring pathways, such as oxidative phosphorylation, and see what genes are most important in gene set enrichment. To confirm these results, metabolic cage data looking into energy usage and expenditure will be extremely insightful. Adding a high fat diet to this system would also be interesting to see how the carb-lipid-protein metabolic homeostasis might shift in the GF mice.
# ggarrange(plot_res$time_all$hm$immune, plot_res$time_all$hm$cancer, plot_res$time_all$hm$ox_phos) %>%
# annotate_figure(top = text_grob("Hallmark", face = "bold", size = 14),
# left = text_grob("pathways", rot = 90),
# bottom = text_grob("Normalized Enrichment Score"))
# f_res$time_1$hm %>% .[rowSums(sapply(., '%in%', setdiff(.$pathway, f_res$time_all$GO$pathway))) > 0,] %>% select(-leadingEdge, -ES, -nMoreExtreme) %>%
# DT::datatable(filter = "bottom")
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggpubr_0.2.5 DT_0.13
## [3] forcats_0.5.0 stringr_1.4.0
## [5] purrr_0.3.3 readr_1.3.1
## [7] tidyr_1.0.2 tibble_2.1.3
## [9] ggplot2_3.3.0 tidyverse_1.3.0
## [11] fgsea_1.10.1 Rcpp_1.0.4.6
## [13] viridis_0.5.1 viridisLite_0.3.0
## [15] RColorBrewer_1.1-2 curl_4.3
## [17] biomaRt_2.40.5 pheatmap_1.0.12
## [19] vsn_3.52.0 GO.db_3.8.2
## [21] org.Mm.eg.db_3.8.2 AnnotationDbi_1.46.1
## [23] AnnotationHub_2.16.1 BiocFileCache_1.8.0
## [25] dbplyr_1.4.2 genefilter_1.66.0
## [27] magrittr_1.5 dplyr_0.8.5
## [29] DESeq2_1.24.0 SummarizedExperiment_1.14.1
## [31] DelayedArray_0.10.0 BiocParallel_1.18.1
## [33] matrixStats_0.56.0 Biobase_2.44.0
## [35] GenomicRanges_1.36.1 GenomeInfoDb_1.20.0
## [37] IRanges_2.18.3 S4Vectors_0.22.1
## [39] BiocGenerics_0.30.0 feather_0.3.5
##
## loaded via a namespace (and not attached):
## [1] snow_0.4-3 readxl_1.3.1
## [3] backports_1.1.5 Hmisc_4.4-0
## [5] fastmatch_1.1-0 splines_3.6.3
## [7] crosstalk_1.1.0.1 digest_0.6.25
## [9] htmltools_0.4.0 fansi_0.4.1
## [11] checkmate_2.0.0 memoise_1.1.0
## [13] cluster_2.1.0 limma_3.40.6
## [15] annotate_1.62.0 modelr_0.1.6
## [17] prettyunits_1.1.1 jpeg_0.1-8.1
## [19] colorspace_1.4-1 blob_1.2.1
## [21] rvest_0.3.5 rappdirs_0.3.1
## [23] haven_2.2.0 xfun_0.12
## [25] crayon_1.3.4 RCurl_1.98-1.1
## [27] jsonlite_1.6.1 survival_3.1-11
## [29] glue_1.3.2 gtable_0.3.0
## [31] zlibbioc_1.30.0 XVector_0.24.0
## [33] scales_1.1.0 DBI_1.1.0
## [35] xtable_1.8-4 progress_1.2.2
## [37] htmlTable_1.13.3 klippy_0.0.0.9500
## [39] foreign_0.8-75 bit_1.1-15.2
## [41] preprocessCore_1.46.0 Formula_1.2-3
## [43] htmlwidgets_1.5.1 httr_1.4.1
## [45] acepack_1.4.1 farver_2.0.3
## [47] pkgconfig_2.0.3 XML_3.99-0.3
## [49] nnet_7.3-13 locfit_1.5-9.2
## [51] labeling_0.3 tidyselect_1.0.0
## [53] rlang_0.4.5 later_1.0.0
## [55] munsell_0.5.0 cellranger_1.1.0
## [57] tools_3.6.3 cli_2.0.2
## [59] generics_0.0.2 RSQLite_2.2.0
## [61] broom_0.5.5 evaluate_0.14
## [63] fastmap_1.0.1 yaml_2.2.1
## [65] knitr_1.28 bit64_0.9-7
## [67] fs_1.3.2 nlme_3.1-144
## [69] mime_0.9 xml2_1.2.5
## [71] compiler_3.6.3 rstudioapi_0.11
## [73] png_0.1-7 interactiveDisplayBase_1.22.0
## [75] ggsignif_0.6.0 affyio_1.54.0
## [77] reprex_0.3.0 geneplotter_1.62.0
## [79] stringi_1.4.6 lattice_0.20-40
## [81] Matrix_1.2-18 vctrs_0.2.4
## [83] pillar_1.4.3 lifecycle_0.2.0
## [85] BiocManager_1.30.10 cowplot_1.0.0
## [87] data.table_1.12.8 bitops_1.0-6
## [89] httpuv_1.5.2 R6_2.4.1
## [91] latticeExtra_0.6-29 affy_1.62.0
## [93] promises_1.1.0 gridExtra_2.3
## [95] assertthat_0.2.1 withr_2.1.2
## [97] GenomeInfoDbData_1.2.1 hms_0.5.3
## [99] grid_3.6.3 rpart_4.1-15
## [101] rmarkdown_2.1 shiny_1.4.0.2
## [103] lubridate_1.7.4 base64enc_0.1-3
smanzoor@uchicago.edu